## ------------------------
## REPLICATION CODE FOR:
## Unemployment Insurance, Risk, and the Acquisition of Specific Skills: An Experimental Approach
## AUTHORS:
## 1. John S. Ahlquist
## 2. Ben Ansell
## ------------------------

# setup -----

rm(list=ls());gc()

## set seed
set.seed(5879521)  

library(here) # 1.0.1
 
##checkpoint dummy file
cat("library(foreign)", "library(MASS)",
    "library(mediation)",
    "library(broom)",
    "library(GGally)",
    "library(ggpubr)",
    "library(plyr)",
    "library(apsrtable)",
    "library(stargazer)",
    "library(readstata13)",
    "library(PropCIs)",
    "library(DescTools)",
    "library(clusterSEs)",
    "library(tidyverse)",
    "library(cobalt)",
    "library(modelsummary)",
    sep="\n",
    file=here("AhlquistAnsell_JPIPE_ReplicationArchive", "checkpoint_dummy.R")) 

## libraries
library(checkpoint) # version 1.0.2
checkpoint("2023-04-01")
library(foreign) # version 0.8.84
library(MASS) # version 7.3.60
library(mediation) # version 4.5.0
library(broom) # version 1.0.4
library(GGally) # version 2.1.2
library(ggpubr) # version 0.6.0
library(stargazer) # version 5.2.3
library(readstata13) # version 0.10.1
library(PropCIs) # version 0.3.0
library(DescTools) # version 0.99.49
library(clusterSEs) # version 2.6.5
library(tidyverse) # version 2.0.0
library(cobalt) # version 4.5.1
library(modelsummary) # version 1.4.2
library(sandwich) # version 3.0.2
library(dataMaid) # version 1.4.1
options("modelsummary_format_numeric_latex" = "plain")


packageVersion("foreign") 

# Figure 1: 2016 employment rate and harmonized unemployment rate against the UI net replacement rate across OECD countries. -------
# The rate (x-axis) is for the initial quarter of an unemployment spell for a married worker who earned the average wage with an out-of-work spouse and two children.
# input: NRRbyEmpOECD.csv
# output: scatterplots of countries and a bivariate best fit line. 

nrr.dat<-read.csv(here("AhlquistAnsell_JPIPE_ReplicationArchive","NRRbyEmpOECD.csv"))

# create Figure 1.a
pdf(here("article", "JPIPE","FinalAcceptedVersion", 
         "figure1a.pdf"))
plot(nrr.dat$NRR_APW_married2c_initial,nrr.dat$employment.rate, col="white", bty="n",
     las=1, xlab="Net UI Replacement Rate (2016)", ylab="Employment Rate (2016)", ylim=c(40,85))
text(nrr.dat$NRR_APW_married2c_initial,nrr.dat$employment.rate, nrr.dat$wbcode)
abline(lm(nrr.dat$employment.rate~nrr.dat$NRR_APW_married2c_initial))
dev.off()

# create Figure 1.b
pdf(here("article", "JPIPE","FinalAcceptedVersion", 
         "figure1b.pdf"))
plot(nrr.dat$NRR_APW_married2c_initial,nrr.dat$HUR, col="white", bty="n", 
     las=1, xlab="Net UI Replacement Rate (2016)", 
     ylab="Harmonized Unemployment Rate (2016)", ylim=c(0,30))
text(nrr.dat$NRR_APW_married2c_initial,nrr.dat$HUR, nrr.dat$wbcode)
abline(lm(nrr.dat$HUR~nrr.dat$NRR_APW_married2c_initial))
dev.off()


# Data Construction -------
# input: AhlquistAnsell_UIskills_cleaned.csv
# output: exp.dat (used for the rest of analysis)

exp.dat<-read.csv(here("AhlquistAnsell_JPIPE_ReplicationArchive", "AhlquistAnsell_UIskills_cleaned.csv"),
                  header=T)
exp.dat$employed<-TRUE
exp.dat$employed[!(exp.dat$emp %in% c("Employed for wages", "Self-employed"))]<-FALSE
exp.dat$college <- FALSE
exp.dat$college[exp.dat$edu %in% c("Bacherlor's degree", "College (Postgraduate Degree)",
                                   "College (Undergraduate Degree)", "Postgraduate degree")]<-TRUE
exp.dat$white<-exp.dat$race=="White"

## Reordering the treatment factor to go from Low to High on both variables in sensible way
exp.dat$treat<-factor(exp.dat$treat, levels = c("LowNone", "LowMinimal", "LowGenerous", "HiNone", "HiMinimal", "HiGenerous"))
summary(exp.dat$treat)
exp.dat$clusters <- as.factor(paste(exp.dat$treat, exp.dat$context, sep=""))

exp.dat$unemprate<-factor(exp.dat$unemprate, levels=c("low", "high"))
exp.dat$UI<-factor(exp.dat$UI, levels = c("none", "minimal", "generous"))
table(exp.dat$context, exp.dat$treat)

## create numerical versions of treatment quantities
exp.dat$UI.num<-0
exp.dat$UI.num[exp.dat$UI=="minimal"]<-0.25
exp.dat$UI.num[exp.dat$UI=="generous"]<-0.75

exp.dat$unemprate.num<-.1
exp.dat$unemprate.num[exp.dat$unemprate=="high"]<-0.25

exp.dat$UI.num<-0
exp.dat$UI.num[exp.dat$UI=="minimal"]<-0.25
exp.dat$UI.num[exp.dat$UI=="generous"]<-0.75

exp.dat$unemprate.num<-.1
exp.dat$unemprate.num[exp.dat$unemprate=="high"]<-0.25

## HH income cleaning
exp.dat$hhinc<-1
exp.dat$hhinc[exp.dat$htotinc == ""] <- NA
exp.dat$hhinc[exp.dat$htotinc %in% c("Â£10-25k", "$15-25k") ] <- 2
exp.dat$hhinc[exp.dat$htotinc %in% c("Â£25-50k", "$25-50k") ] <- 3
exp.dat$hhinc[exp.dat$htotinc %in% c("Â£50-100k", "$50-100k") ] <- 4
exp.dat$hhinc[grepl("More", exp.dat$htotinc)] <- 5

## unemployed IRL
exp.dat$unemployed.irl<-0
exp.dat$unemployed.irl[exp.dat$emp %in% c("Out of work and looking for work", "Other")] <-1

## CRRA parameter implied from risk preference elicitation
### gamble calculations
crra<-function(x,theta){x^(1-theta)/(1-theta)}
p<-.5
####expected values
E.g1<-80
E.g2<-p*120 + (1-p)*60
E.g3<-p*160 + (1-p)*40
E.g4<-p*200 - (1-p)*20
E.g5<-p*240 - (1-p)*0
u.g1<-function(theta){p*crra(80,theta) + (1-p)*crra(80,theta)} 
u.g2<-function(theta){p*crra(120,theta) + (1-p)*crra(60,theta)} 
u.g3<-function(theta){p*crra(160,theta) + (1-p)*crra(40,theta)} 
u.g4<-function(theta){p*crra(200,theta) + (1-p)*crra(20,theta)} 
u.g5<-function(theta){p*crra(240,theta) + (1-p)*crra(0,theta)}

#### finding risk aversion parameter that makes you indifferent between 2 gambles
theta.1<-uniroot(function(theta)u.g1(theta)-u.g2(theta), c(0,3))$root #2 
theta.2<-uniroot(function(theta)u.g2(theta)-u.g3(theta), c(0,3))$root #0.67
theta.3<-uniroot(function(theta)u.g3(theta)-u.g4(theta), c(0,3))$root # .38
theta.4<-uniroot(function(theta)u.g4(theta)-u.g5(theta), c(0,3))$root #0.20

exp.dat$risk.param<-NA
exp.dat$risk.param[grepl("80 ECU", exp.dat$risk)]<-theta.1
exp.dat$risk.param[grepl("120 ECU", exp.dat$risk)]<-mean(c(theta.1, theta.2))
exp.dat$risk.param[grepl("160 ECU", exp.dat$risk)]<-mean(c(theta.2, theta.3))
exp.dat$risk.param[grepl("200 ECU", exp.dat$risk)]<-mean(c(theta.3, theta.4))
exp.dat$risk.param[grepl("240 ECU", exp.dat$risk)]<-theta.4

exp.dat$risk.num<-NA
exp.dat$risk.num[grepl("80 ECU", exp.dat$risk)]<-4
exp.dat$risk.num[grepl("120 ECU", exp.dat$risk)]<-3
exp.dat$risk.num[grepl("160 ECU", exp.dat$risk)]<-2
exp.dat$risk.num[grepl("200 ECU", exp.dat$risk)]<-1
exp.dat$risk.num[grepl("240 ECU", exp.dat$risk)]<-0

# create codebook of exp.dat (cleaned dataset) -----
# makeDataReport(exp.dat, file = "AhlquistAnsell_JPIPE_ReplicationArchive/codebook.Rmd", replace = TRUE)  # Produce data report by DataMaid

# Figure 5: Standardized differences in means for covariates -----
# Each point represents the maximum difference between the mean of a covariate in one treatment condition and the mean of that covariate in all other experimental conditions.

baltab<-bal.tab(treat ~ context + sex + age + white + college + employed + hhinc + risk.param, #+ duration, 
                data=exp.dat, s.d.denom="pooled", pairwise=FALSE)

pdf(here("article", "JPIPE","FinalAcceptedVersion", 
         "figure5.pdf"))
plot(baltab, stars = "raw")
dev.off()   

# subset the dataset to separate `context` (lab or mTurk)
exp.dat.lab<-subset(exp.dat, context=="lab")
exp.dat.mturk<-subset(exp.dat, context=="mTurk")

# Figure 4: The proportion of subjects giving up 1000ECU to acquire slider-specific skill, as a function of experimental treatments. ------
## Figure 4.b: pooled across UI levels
i_tab.UI<-table(exp.dat$investment, exp.dat$UI)
tab.UI<-margin.table(table(exp.dat$investment, exp.dat$UI),2)
sum.treat.UI<-data.frame(BinomCI(i_tab.UI[2,], tab.UI))
xs<-c(0,.25,.75)

pdf(file = here(
  "article", "JPIPE","FinalAcceptedVersion", 
  "figure4b.pdf"),
  width=6.0, height=5.3, colormodel = "cmyk")
plot(xs, (sum.treat.UI$est[1:3]), type="b", las=1,
     xlim=c(0,1), ylim=c(.3,.9), bty="n", ylab="Proportion Making Investment",
     xlab="Generosity", xaxt="n", lwd=3, cex=1.5, col="#14284B")
axis(1, at=c(0,.25,.75))
segments(xs,(sum.treat.UI$lwr.ci[1:3]),xs,(sum.treat.UI$upr.ci[1:3]),
         col="#14284B") 
dev.off()
rm(xs)

## Figure 4.a: pooled across unemployment risk
i_tab.unemp<-table(exp.dat$investment, exp.dat$unemprate)
tab.unemp<-margin.table(table(exp.dat$investment, exp.dat$unemprate),2)
sum.treat.unemp<-data.frame(BinomCI(i_tab.unemp[2,], tab.unemp))
xs<-c(0.1, 0.25)
pdf(file= here(
  "article", "JPIPE","FinalAcceptedVersion", 
  "figure4a.pdf"),
  width=6.0, height=5.3, colormodel = "cmyk")
plot(xs, (sum.treat.unemp$est[1:2]), type="b", las=1,
     xlim=c(0.1,.3), ylim=c(.5,.75), bty="n", ylab="Proportion Making Investment",
     xlab="Unemployment Risk", xaxt="n", lwd=3, cex=1.5, col="#14284B")
axis(1, at=c(.1,.25), labels = c("10%", "25%"))
segments(xs,(sum.treat.unemp$lwr.ci[1:2]),xs,(sum.treat.unemp$upr.ci[1:2]),
         col="#14284B") 
dev.off()
rm(xs)

## Figure 4.c: Investment decision
invest.tab<-table(exp.dat$investment, exp.dat$treat)
tab.n<-margin.table(table(exp.dat$investment, exp.dat$treat),2)
sum.treat<-data.frame(BinomCI(invest.tab[2,], tab.n))
xs<-c(0,.25,.75)

pdf(file= here(
  "article", "JPIPE","FinalAcceptedVersion", 
  "figure4c.pdf"),
  width=6, height=5.3, colormodel = "cmyk")
plot(xs, (sum.treat$est[4:6]), type="b", las=1,
     xlim=c(0,1), ylim=c(.3,.9), bty="n", ylab="proportion",
     xlab="Generosity", xaxt="n", lwd=3, cex=1.5, col="#14284B")
axis(1, at=c(0,.25,.75))
lines(xs+.02, (sum.treat$est[1:3]), type="b", col="darkorange2", pch=2,
      cex=1.5, lwd=3)
segments(xs,(sum.treat$lwr.ci[4:6]),xs,(sum.treat$upr.ci[4:6]),
         col="#14284B")
segments(xs+.02,(sum.treat$lwr.ci[1:3]),xs+.02,(sum.treat$upr.ci[1:3]),
         col="darkorange2")

text(x=.3, y=.8,labels="25% unemployment", col="#14284B")
text(x=.3, y=.45,labels="10% unemployment", col="darkorange2")
dev.off()
rm(xs)

# Table 4: p-values for differences in the proportion of subjects choosing to invest, adjusted to control the false discovery rate for multiple comparisons.----

## Using BH correction to control fdr for multiple comparisons
investpair.bh<-pairwise.prop.test(invest.tab[2,], tab.n,"BH")
investpair<-pairwise.prop.test(invest.tab[2,], tab.n,"none")

p.tab<-round(investpair$p.value,3)
p.tab.bh<-round(investpair.bh$p.value,3)

for(i in 1:nrow(p.tab)){
  for(j in 1:ncol(p.tab)){
    p.tab[i,j]<-paste(p.tab.bh[i,j],"(",p.tab[i,j],")", sep="")
  }}
p.tab[p.tab=="NA(NA)"]<-""
stargazer(p.tab, style="apsr", digits =2,
          title = "$p$-values for differences in the proportion of subjects choosing to invest, adjusted to control the false discovery rate for multiple comparisons.  Values in parentheses are uncorrected $p$-values ",
          out = here(
            "article", "JPIPE","FinalAcceptedVersion", 
            "table4.tex"),
          label = "tab:investpvals", font.size = "footnotesize")

# investment regression models -----
## logit
inv.1<-glm(investment=="yes" ~ treat, family=binomial, data=exp.dat)

## evidence of linear effect of UI, so entering linearly
inv.1n<-glm(investment=="yes" ~ unemprate.num + UI.num, 
            family=binomial, data=exp.dat)

inv.1ni<-glm(investment=="yes" ~ unemprate.num*UI.num, 
             family=binomial, data=exp.dat)

inv.1ci<-glm(investment=="yes" ~ unemprate.num*context + UI.num*context, 
             family=binomial, data=exp.dat)

inv.2<-update(inv.1, .~.+ context + I(sex=="Female") + age + I(race=="White") +
                I(area_raised=="Urban") + college+ hhinc + 
                risk.param + unemployed.irl)

inv.2n<-update(inv.1n, .~.+ context + I(sex=="Female") + age + I(race=="White") +
                 I(area_raised=="Urban") + college + hhinc + 
                 risk.param + unemployed.irl)

inv.2n.raint<-update(inv.2n, .~.+ risk.param*unemprate.num) #testing for interaction between risk and risk aversion.

### OLS/LPM
inv.lm.1<-glm(investment=="yes" ~ treat, data=exp.dat)
inv.lm.1n<-glm(investment=="yes" ~ unemprate.num + UI.num, data=exp.dat)
inv.lm.1ni<-glm(investment=="yes" ~ unemprate.num*UI.num, data=exp.dat)
inv.lm.1ci<-glm(investment=="yes" ~ unemprate.num*context + 
                  UI.num*context, data=exp.dat)
inv.lm.2n<-update(inv.lm.1n, .~. + context + I(sex=="Female") + age + 
                    I(race=="White") + I(area_raised=="Urban") + 
                    college + hhinc + risk.param + unemployed.irl)

coef_renames <- c("treatLowMinimal" = "low-minimal",
                  "treatLowGenerous" = "low-generous",
                  "treatHiNone" = "high-none",
                  "treatHiMinimal" = "high-minimal",
                  "treatHiGenerous" = "high-generous",
                  "unemprate.num"= "unemployment rate",
                  "UI.num" = "UI generousity",
                  "contextmTurk" = "mTurk",
                  "unemprate.num:UI.num" = "generosity $\times$ rate",
                  "unemprate.num:contextmTurk" = "rate $\times$ mTurk",
                  "contextmTurk:UI.num" = "generosity $\times$ mTurk")

coef_renames_full <- c("treatLowMinimal" = "low-minimal",
                       "treatLowGenerous" = "low-generous",
                       "treatHiNone" = "high-none",
                       "treatHiMinimal" = "high-minimal",
                       "treatHiGenerous" = "high-generous",
                       "unemprate.num"= "unemployment rate",
                       "UI.num" = "UI generousity",
                       "contextmTurk" = "mTurk",
                       "unemprate.num:UI.num" = "generosity $\times$ rate",
                       "unemprate.num:contextmTurk" = "rate $\times$ mTurk",
                       "contextmTurk:UI.num" = "generosity $\times$ mTurk",
                       'I(sex == "Female")TRUE' = "female",
                       "age" = "age",
                       'I(race == "White")TRUE' = "white",
                       'I(area_raised == "Urban")TRUE' = "urban",
                       "collegeTRUE" = "college degree",
                       "hhinc" = "income",
                       "risk.param" = "risk aversion",
                       "unemployed.irl" = "unemployed"
)


# Table 7: Logistic Regression for Investment Choice ------
modelsummary(list("base (categorical)" = inv.1, 
                  "base (linear)" = inv.1n,
                  "interacted" = inv.1ni,
                  "context" = inv.1ci, 
                  "covariates" = inv.2n),
             output = here(
               "article", "JPIPE","FinalAcceptedVersion",
               "table7.tex"),
             title = "Logistic Regression for Investment Choice \\label{tab:invest_main_logit}",
             fmt=2, stars = TRUE, coef_map = coef_renames,
             # Specify how to robustify/cluster the model
             vcov = function(x) vcovCL(x, cluster = ~ clusters),
             # Get rid of other coefficients and goodness-of-fit (gof) stats
             gof_omit = 'DF|Deviance|R2|AIC|RMSE|logLik|se_type'
)

# Table 3: Linear Regression for Investment Choice -----
modelsummary(list("base (categorical)" = inv.lm.1, 
                  "base (linear)" = inv.lm.1n,
                  "interacted" = inv.lm.1ni, 
                  "context" = inv.lm.1ci, 
                  "covariates" = inv.lm.2n),
             output = here(
               "article", "JPIPE","FinalAcceptedVersion",
               "table3.tex"),
             title = "Linear Regression for Investment Choice \\label{tab:invest_main_lm}",
             fmt=2, stars = TRUE, coef_map = coef_renames,
             # Specify how to robustify/cluster the model
             vcov = function(x) vcovCL(x, cluster = ~ clusters),
             # Get rid of other coefficients and goodness-of-fit (gof) stats
             gof_omit = 'DF|Deviance|AIC|RMSE|logLik|se_type'
)

# Table 5: Linear Regression for Investment Choice (detailed results) -----
modelsummary(list("base (categorical)" = inv.lm.1, 
                  "base (linear)" = inv.lm.1n,
                  "interacted" = inv.lm.1ni, 
                  "context" = inv.lm.1ci, 
                  "covariates" = inv.lm.2n),
             output = here(
               "article", "JPIPE","FinalAcceptedVersion",
               "table5.tex"),
             title = "Linear Regression for Investment Choice (detailed results) \\label{tab:invest_main_lm_cov}",
             fmt=2, stars = TRUE, coef_map = coef_renames_full,
             # Specify how to robustify/cluster the model
             vcov = function(x) vcovCL(x, cluster = ~ clusters),
             # Get rid of other coefficients and goodness-of-fit (gof) stats
             gof_omit = 'DF|Deviance|AIC|RMSE|logLik|se_type'
)


# classical SEs (for appendix) ----
## Table 8: Logistic Regression for Investment Choice (classical standard errors) -----
modelsummary(list("base (categorical)" = inv.1, 
                  "base (linear)" = inv.1n,
                  "interacted" = inv.1ni,
                  "context" = inv.1ci, 
                  "covariates" = inv.2n),
             output = here(
               "article", "JPIPE","FinalAcceptedVersion",
               "table8.tex"),
             title = "Logistic Regression for Investment Choice (classical standard errors) \\label{tab:invest_logit_classicSE}",
             fmt=2, stars = TRUE, coef_map = coef_renames,
             # Get rid of other coefficients and goodness-of-fit (gof) stats
             gof_omit = 'DF|Deviance|R2|AIC|RMSE|logLik|se_type'
)

## table 6: Linear Regression for Investment Choice (robust standard errors) -----
modelsummary(list("base (categorical)" = inv.lm.1, 
                  "base (linear)" = inv.lm.1n,
                  "interacted" = inv.lm.1ni,
                  "context" = inv.lm.1ci, 
                  "covariates" = inv.lm.2n),
             output = here(
               "article", "JPIPE","FinalAcceptedVersion",
               "table6.tex"),
             title = "Linear Regression for Investment Choice (robust standard errors) \\label{tab:invest_lm_classicSEs}",
             fmt=2, stars = TRUE, coef_map = coef_renames,
             # Specify how to robustify/cluster the model
             vcov = "robust",
             # Get rid of other coefficients and goodness-of-fit (gof) stats
             gof_omit = 'DF|Deviance|AIC|RMSE|logLik|se_type'
)

# Breaking out Investment results by Lab v mTurk for Appendix ----
## LAB SESSIONS ----

invest.tab.lab<-table(exp.dat.lab$investment, exp.dat.lab$treat)
tab.n.lab<-margin.table(table(exp.dat.lab$investment, exp.dat.lab$treat),2)
sum.treat.lab<-data.frame(BinomCI(invest.tab.lab[2,], tab.n.lab))
xs<-c(0,.25,.75)

### Figure 6.a of basic investment results -----
pdf(
  file = here("article", "JPIPE","FinalAcceptedVersion",
              "figure6a.pdf"),
  width=6.8, height=5.45, colormodel = "cmyk")
plot(xs, (sum.treat.lab$est[4:6]), type="b", las=1,
     xlim=c(0,1), ylim=c(.3,.9), bty="n", ylab="proportion",
     xlab="generosity", xaxt="n", lwd=3, cex=1.5, col="#14284B", main= "Lab subjects")
axis(1, at=c(0,.25,.75))
lines(xs+.02, (sum.treat.lab$est[1:3]), type="b", col="darkorange2", pch=2,
      cex=1.5, lwd=3)
segments(xs,(sum.treat.lab$lwr.ci[4:6]),xs,(sum.treat.lab$upr.ci[4:6]),
         col="#14284B")
segments(xs+.02,(sum.treat.lab$lwr.ci[1:3]),xs+.02,(sum.treat.lab$upr.ci[1:3]),
         col="darkorange2")
text(x=.5, y=.8,labels="25% unemployment", col="#14284B")
text(x=.5, y=.65,labels="10% unemployment", col="darkorange2")
dev.off()

### Table 9: p-values for differences in the proportion of subjects choosing to invest, lab subjects only. -----
investpair.bh.lab<-pairwise.prop.test(invest.tab.lab[2,], tab.n.lab,"BH")
investpair.lab<-pairwise.prop.test(invest.tab.lab[2,], tab.n.lab,"none")
p.tab.lab<-round(investpair.lab$p.value,2)
p.tab.bh.lab<-round(investpair.bh.lab$p.value,2)

for(i in 1:nrow(p.tab.lab)){
  for(j in 1:ncol(p.tab.lab)){
    p.tab.lab[i,j]<-paste(p.tab.bh.lab[i,j],"(",p.tab.lab[i,j],")", sep="")
  }}
p.tab.lab[p.tab.lab=="NA(NA)"]<-""
stargazer(p.tab.lab, style="apsr", 
          title = "$p$-values for differences in the proportion of subjects choosing to invest, lab subjects only.  Values adjusted using the Bejamini-Hotchberg correction to control false discovery rate for multiple comparisons.  Values in parentheses are unadjusted $p$-values ",
          out = here("article", "JPIPE","FinalAcceptedVersion", 
                     "table9.tex"),
          label = "tab:investpvals_lab")

## MTURK SESSIONS -----

invest.tab.mturk<-table(exp.dat.mturk$investment, exp.dat.mturk$treat)
tab.n.mturk<-margin.table(table(exp.dat.mturk$investment, exp.dat.mturk$treat),2)
sum.treat.mturk<-data.frame(BinomCI(invest.tab.mturk[2,], tab.n.mturk))
xs<-c(0,.25,.75)

### Figure 6.b of basic investment results ----
pdf(
  file = here("article", "JPIPE","FinalAcceptedVersion", 
              "figure6b.pdf"),
  width=6.8, height=5.45, colormodel = "cmyk")
plot(xs, (sum.treat.mturk$est[4:6]), type="b", las=1,
     xlim=c(0,1), ylim=c(.3,.9), bty="n", ylab="proportion",
     xlab="generosity", xaxt="n", lwd=3, cex=1.5, col="#14284B", main= "mTurk subjects")
axis(1, at=c(0,.25,.75))
lines(xs+.02, (sum.treat.mturk$est[1:3]), type="b", col="darkorange2", pch=2,
      cex=1.5, lwd=3)
segments(xs,(sum.treat.mturk$lwr.ci[4:6]),xs,(sum.treat.mturk$upr.ci[4:6]),
         col="#14284B")
segments(xs+.02,(sum.treat.mturk$lwr.ci[1:3]),xs+.02,(sum.treat.mturk$upr.ci[1:3]),
         col="darkorange2")
text(x=.5, y=.75,labels="25% unemployment", col="#14284B")
text(x=.5, y=.5,labels="10% unemployment", col="darkorange2")
dev.off()

### Table 10: Pairwise comparisons MTURK -----
investpair.bh.mturk<-pairwise.prop.test(invest.tab.mturk[2,], tab.n.mturk,"BH")
investpair.mturk<-pairwise.prop.test(invest.tab.mturk[2,], tab.n.mturk,"none")


p.tab.mturk<-round(investpair.mturk$p.value,2)
p.tab.bh.mturk<-round(investpair.bh.mturk$p.value,2)


for(i in 1:nrow(p.tab.mturk)){
  for(j in 1:ncol(p.tab.mturk)){
    p.tab.mturk[i,j]<-paste(p.tab.bh.mturk[i,j],"(",p.tab.mturk[i,j],")", sep="")
  }}
p.tab.mturk[p.tab.mturk=="NA(NA)"]<-""
stargazer(p.tab.mturk, style="apsr", 
          title = "$p$-values for differences in the proportion of subjects choosing to invest, mTurk subjects only.  Values adjusted to control the false discovery rate for multiple comparisons.  Values in parentheses are uncorrected $p$-values ",
          out = here("article", "JPIPE","FinalAcceptedVersion", 
                     "table10.tex"),
          label = "tab:investpvals_mturk")

### Table 11: Regression table, separating by context -----
lab.lm.1n<-update(inv.lm.1n, data = exp.dat.lab)
lab.lm.2n<-update(inv.lm.2n, .~. - context, data = exp.dat.lab)
mturk.lm.1n<-update(inv.lm.1n, data = exp.dat.mturk)
mturk.lm.2n<-update(inv.lm.2n, .~. - context, data = exp.dat.mturk)

samp <- tribble(~term, ~base, ~covariates, ~base, ~covariates,
                "Context", "lab", "lab", "mTurk","mTurk")
coef_bycontext<-coef_renames[c(6,7)]
modelsummary(list(
  "base" = lab.lm.1n, 
  "covariates" = lab.lm.2n,
  "base" = mturk.lm.1n, 
  "covariates" = mturk.lm.2n),
  output = here(
    "article", "JPIPE","FinalAcceptedVersion",
    "table11.tex"),
  title = "Linear Regression for Investment Choice by experimental context \\label{tab:tab:investContext_num_lm}",
  fmt=2, stars = TRUE, coef_map = coef_bycontext,
  # Specify how to robustify/cluster the model
  vcov = "robust",
  # Get rid of other coefficients and goodness-of-fit (gof) stats
  gof_omit = 'DF|Deviance|AIC|RMSE|logLik|se_type',
  add_rows = samp
)

# Waiting Time ------------------------------------------------------------

wait.dat<-subset(exp.dat, unemp.exp==TRUE) #only subjects experiencing unemployment

wait2.tab<-table(wait.dat$waiting2, wait.dat$treat)
tab.n<-margin.table(table(wait.dat$waiting2, wait.dat$treat),2)
sum.treat<-data.frame(BinomCI(wait2.tab[2,], tab.n))
xs<-c(0,.25,.75)

## Figure 8: Proportion of lab respondents having had at least one spell of unemployment who chose to extend their unemployment by rejecting a round of work at least once. -----
pdf(file= here(
  "article", "JPIPE","FinalAcceptedVersion",
  "figure8.pdf"),
  width=6, height=5.3, colormodel = "cmyk")

plot(xs, (sum.treat$est[4:6]), type="b", las=1,
     xlim=c(0,1), ylim=c(0,.5), bty="n", ylab="proportion",
     xlab="Generosity", xaxt="n", lwd=3, cex=1.5, col="#14284B")
axis(1, at=c(0,.25,.75))
lines(xs+.02, (sum.treat$est[1:3]), type="b", col="darkorange2", pch=2,
      cex=1.5, lwd=3)
segments(xs,(sum.treat$lwr.ci[4:6]),xs,(sum.treat$upr.ci[4:6]),
         col="#14284B")
segments(xs+.02,(sum.treat$lwr.ci[1:3]),xs+.02,(sum.treat$upr.ci[1:3]),
         col="darkorange2")

text(x=.6, y=.4,labels="25% unemployment", col="#14284B")
text(x=.42, y=.1,labels="10% unemployment", col="darkorange2")
dev.off()

# Figure 7: The distribution of the number of foregone employment offers by unemployment risk among lab subjects who have at least one spell of unemployment. -----
pdf( file= here(
  "article", "JPIPE","FinalAcceptedVersion",
  "figure7.pdf"),
  width=6, height=5.3, colormodel = "cmyk")

hlue<-hist(wait.dat$tot.rejemp[wait.dat$unemprate=="low"], breaks= 0:6, plot=F, right=F)
hhue<-hist(wait.dat$tot.rejemp[wait.dat$unemprate=="high"], breaks= 0:6, plot=F, right=F)
hgui<-hist(wait.dat$tot.rejemp[wait.dat$UI=="generous"], breaks= 0:5, plot=F)
hmui<-hist(wait.dat$tot.rejemp[wait.dat$UI=="minimal"], breaks= 0:5, plot=F)
hnui<-hist(wait.dat$tot.rejemp[wait.dat$UI=="none"], breaks= 0:5, plot=F)

plot.df<-data.frame(
  proportion = c(hhue$density,hlue$density),
  unemployment = c(rep("high",length(hhue$density)), 
                   rep("low",length(hlue$density))
  ),
  waiting = c(0:(max(hlue$breaks)-1),0:(max(hlue$breaks)-1))
)

ggplot(plot.df, aes(waiting, proportion)) +
  geom_linerange(
    aes(x = waiting, ymin = 0, ymax = proportion, group = unemployment), 
    color = "lightgray", size = 1.5,
    position = position_dodge(0.3)
  ) +
  geom_point(
    aes(color = unemployment),
    position = position_dodge(0.3), size = 3
  ) +
  scale_color_manual(values = c("#0073C2FF", "#EFC000FF"))+
  theme_pubclean()
dev.off()


# Table 12: Regression Analyses for Waiting -----

wait1.pois<-glm(tot.rejemp ~ UI + unemprate + log(offers), 
                data = subset(wait.dat, offers>0), family = "poisson")

wait1.pois.cov<-update(wait1.pois, .~. + I(sex=="Female") + age + 
                         I(race=="White") + I(area_raised=="Urban") + college + risk.param +
                         unemployed.irl + hhinc) 

wait2.reg<-glm(waiting2~UI + unemprate , family=binomial, data= wait.dat)
wait2a.reg<-glm(waiting2~UI + unemprate + offers, family=binomial, data= wait.dat)
wait2a.reg.cov<-update(wait2a.reg, .~. +I(sex=="Female") + age + I(race=="White") +
                         I(area_raised=="Urban") + college + risk.param +
                         unemployed.irl + hhinc)

stargazer(wait2.reg, wait2a.reg, wait2a.reg.cov, wait1.pois, wait1.pois.cov, 
          style="apsr", out = here(
            "article", "JPIPE","FinalAcceptedVersion",
            "table12.tex"),
          title = "Regression Analyses for Waiting", digits = 2,
          label = "tab:waiting", dep.var.caption = "", 
          dep.var.labels.include = FALSE, dep.var.labels = c("", ""),
          column.labels = c("Ever Rejected", "Number of Rejected "), 
          column.separate =c(3,2),
          covariate.labels=c("Minimal UI",  "Generous UI", "High risk", 
                             "emp. offers", "Female", "Age", "White", 
                             "Urban", "College", "Risk aversion", 
                             "unemployed", "income", "log cont", "Constant"),
          keep.stat=c("n", "aic"), font.size = "footnotesize")

# Figure 12: Mediation analysis of skill investment on waiting behavior ----

wait.dat$inv.num<-as.numeric(as.factor(wait.dat$investment))-1
out.fit<-glm(waiting2~ inv.num + unemprate.num +UI.num + offers, 
             data=wait.dat, family = binomial)
obs<-row.names(out.fit$model)
med.fit<-glm(inv.num~ unemprate.num + UI.num +college+risk.param,
             family=binomial, 
             data=wait.dat[row.names(wait.dat)%in% obs,])
med.out.risk<-mediate(med.fit, out.fit, treat="unemprate.num", 
                      mediator="inv.num", sims=1000)
med.out.UI<-mediate(med.fit, out.fit, treat="UI.num", 
                    mediator="inv.num", sims=1000)
pdf(file = here(
  "article", "JPIPE","FinalAcceptedVersion",
  "figure12.pdf"),
  width=9, height=5, colormodel = "cmyk")
par(mfrow=c(1,2))
plot(med.out.risk, treatment = "treated", bty="n", lwd=2,
     main="Effects of Unemployment Risk on Waiting: \n Mediated by Investment")
plot(med.out.UI, treatment = "treated", bty="n", lwd=2,
     main="Effects of Generosity on Waiting: \n Mediated by Investment")
dev.off()

# Task Effort Appendix F -----------------------------------------------------

# Box plots 

# Figure 9: Effort by Unemployment Risk -----
p <- ggboxplot(exp.dat, x = "unemprate", y = "effort2", 
               color="unemprate",  add="jitter")

p + stat_compare_means(method = "t.test", hjust=-0.7)+
  xlab("Unemployment Rate")+
  ylab("Average Effort per Round")+
  theme(legend.position="none")

ggsave(here(
  "article", "JPIPE","FinalAcceptedVersion",
  "figure9.pdf"),
  width = 6,
  height = 5)

# Figure 10: Effort by Unemployment Insurance ----
p <- ggboxplot(exp.dat, x = "UI", y = "effort2", color="UI",  add="jitter")

p + scale_x_discrete(breaks=c("generous","minimal","none"),
                     labels=c("Generous", "Minimal", "None"))+
  xlab("Unemployment Insurance")+
  ylab("Average Effort per Round")+
  theme(legend.position="none")

ggsave(here(
  "article", "JPIPE","FinalAcceptedVersion",
  "figure10.pdf"),
  width = 6,
  height = 5)

# Figure 11: Effort by Investment Choice -----
p <- ggboxplot(exp.dat, x = "investment", y = "effort2", 
               color="investment",  add="jitter")

p +stat_compare_means(method = "t.test", hjust=-0.7)+
  xlab("Made Investment")+
  ylab("Average Effort per Round")+
  theme(legend.position="none")

ggsave(here(
  "article", "JPIPE","FinalAcceptedVersion",
  "figure11.pdf"),
  width = 6,
  height = 5)

## Table 13: OLS regressions on average task effort ----

eff.1<-lm(effort2 ~ treat, data=exp.dat)
eff.1a<-lm(effort2 ~ treat+investment, data=exp.dat)
eff.1b<-lm(effort2 ~ treat*investment, data=exp.dat)
eff.2<-update(eff.1a, .~.+ context + I(sex=="Female") + age + I(race=="White") +
                area_raised + college + risk.param + hhinc + unemployed.irl)
eff.3<-update(eff.1a, .~.+ investment+ context + I(sex=="Female") + age + I(race=="White") +
                area_raised + college + risk.param + hhinc + unemployed.irl)

efreg1<-lm(effort2 ~ UI.num+unemprate.num, data=exp.dat)

efreg3<-update(efreg1, .~.+ context + I(sex=="Female") + age + I(race=="White") +
                 I(area_raised=="Urban") + college + risk.param + hhinc + unemployed.irl)

stargazer(efreg1, efreg3, style="apsr", type="latex",
          out = here(
            "article", "JPIPE","FinalAcceptedVersion",
            "table13.tex"),
          digits = 2,
          title = "OLS regression on average task effort",
          label = "tab:effort", dep.var.caption = "", dep.var.labels.include = FALSE, 
          dep.var.labels = c("", ""), covariate.labels=c(
            "UI generosity", "Unemployment risk", "mTurk",  
            "Female", "Age", "White", "Urban", "College", "Risk aversion", 
            "income", "unemployed", "Constant"),
          keep.stat=c("n", "adj.rsq"), font.size = "footnotesize")



## Figure 13: Mediation analysis of skill investment on task effort.------

exp.dat$inv.num<-as.numeric(exp.dat$investment)-1
out.fit <- lm(effort2~ inv.num + unemprate.num +UI.num + context + age + risk.param, data=exp.dat)
obs<-row.names(out.fit$model)
med.fit<-glm(inv.num~ unemprate.num + UI.num +college+risk.param, family=binomial, 
             data=exp.dat[row.names(exp.dat)%in% obs,])
med.out.risk<-mediate(med.fit, out.fit, treat="unemprate.num", mediator="inv.num", sims=1000)
med.out.UI<-mediate(med.fit, out.fit, treat="UI.num", mediator="inv.num", sims=1000)
pdf(file= here( "article", "JPIPE","FinalAcceptedVersion",
  "figure13.pdf"),
    width=9, height=5, colormodel = "cmyk")
par(mfrow=c(1,2))
plot(med.out.risk, treatment = "treated", bty="n", lwd=2,
     main="Effects of Unemployment Risk on Effort: \n Mediated by Investment")
plot(med.out.UI, treatment = "treated", bty="n", lwd=2,
     main="Effects of Generosity on Effort: \n Mediated by Investment")
dev.off()

